# matlab代码是最初代码
""""
close all
clear
clc
NUM = 64; 
l = [0, 1]; 
p = [0, 1]; 
n = 1000;   
threshold = 0.95;    


% 生成he 光学数据
LG = LG_modes(NUM, l, p);
datas = generate_random_superpositions(LG, n);

% PCA和绘图
[PC_exp, components_exp, scores_exp, n_min] = PCA(datas, threshold);
plot_results(PC_exp, components_exp, scores_exp, n_min, NUM);


fprintf('\n累计解释方差:\n');
for i = 1:n_min
    fprintf('PC%d: %.1f%%\n', i, components_exp(i));
end
fprintf('前%d个PC累计: %.1f%%\n', n_min, sum(components_exp(1:n_min)));


%% PCA函数
function [PC, components, scores, n_min] = PCA(data, threshold)
    [rows, cols, num_samples] = size(data);
    
    % 重塑数据矩阵
    data_matrix = reshape(data, [rows*cols, num_samples])';
    
    % 中心化数据
    mean_data = mean(data_matrix, 1);
    data_centered = data_matrix - mean_data;
    
    % 计算协方差矩阵
    covariance_matrix = (data_centered' * data_centered) / (num_samples - 1);
    
    % 特征值分解
    [eigenvectors, eigenvalues] = eig(covariance_matrix);
    eigenvalues = diag(eigenvalues);
    
    % 按特征值降序排序
    [eigenvalues, sort_idx] = sort(eigenvalues, 'descend');
    eigenvectors = eigenvectors(:, sort_idx);
    
    % 计算累计解释方差
    total_variance = sum(eigenvalues);
    cumulative_variance = cumsum(eigenvalues) / total_variance;
    
    % 找到达到方差阈值所需的主成分数量
    n_min = find(cumulative_variance >= threshold, 1);
    if isempty(n_min)
        n_min = length(eigenvalues);
    end
    
    % 提取所需的主成分
    PC = eigenvectors(:, 1:n_min);
    scores = data_centered * PC;
    components = 100 * eigenvalues(1:n_min) / total_variance;
end

%% 绘图函数
function plot_results(PC, components, scores, n_min, N)

    tiledlayout('flow','TileSpacing','tight','Padding','tight')
    for i = 1:n_min
        nexttile;
        pc_img = reshape(PC(:, i), [N, N]);
        imagesc(pc_img);
        colormap(jet);
        axis square off;
        title(sprintf('PC%d (%.1f%%)', i, components(i)), 'FontSize', 10);
    end
    
    % 散点图
    nexttile;
    scatter(scores(:, 1), scores(:, 2), 20, 'filled', ...
        'MarkerFaceAlpha', 0.6, 'MarkerEdgeAlpha', 0.6, ...
        'MarkerFaceColor', [0.2, 0.4, 0.8]);
    xlabel('PC1');
    ylabel('PC2');
    grid on;
    axis square
    
    % 方差条形图
    nexttile;
    bar(components(1:min(10, n_min)));
    xlabel('主成分');
    ylabel('解释方差 (%)');
    title('各主成分方差');
    axis square
    
    % 累计方差图
    nexttile;
    cumulative = cumsum(components);
    plot(cumulative, 'o-', 'LineWidth', 2);
    hold on;
    plot([1, n_min], [95, 95], 'r--', 'LineWidth', 1);
    xlabel('主成分数量');
    ylabel('累计解释方差 (%)');
    title('累计解释方差');
    axis square
    ylim([0,100])
end



%% 光学函数
function LG_modes = LG_modes(N, l_values, p_values)
    x = linspace(-3, 3, N);
    [X, Y] = meshgrid(x);
    R = sqrt(X.^2 + Y.^2);
    Phi = atan2(Y, X);
    
    w0 = 1;
    LG_modes = [];
    
    for l = l_values
        for p = p_values
            L = laguerreL_iterative(p, abs(l), 2*R.^2/w0^2);
            amplitude = (sqrt(2)*R/w0).^abs(l) .* L .* exp(-R.^2/w0^2);
            phase = exp(1i*l*Phi);
            mode = amplitude .* phase;
            mode = mode / sqrt(sum(abs(mode(:)).^2));
            LG_modes = cat(3, LG_modes, mode);
        end
    end
end

function result = laguerreL_iterative(p, l, x)
    if p == 0
        result = ones(size(x));
    elseif p == 1
        result = 1 + l - x;
    else
        L_prev2 = ones(size(x));
        L_prev1 = 1 + l - x;
        
        for n = 2:p
            L_current = ((2*n + l - 1 - x) .* L_prev1 - (n + l - 1) * L_prev2) / n;
            L_prev2 = L_prev1;
            L_prev1 = L_current;
        end
        result = L_prev1;
    end
end

function superpositions = generate_random_superpositions(LG_modes, num_samples)
    
    [rows, cols, num_modes] = size(LG_modes);
    superpositions = zeros(rows, cols, num_samples);
    coefficients = zeros(num_samples, num_modes);
    total_intensity=1;
    
    for i = 1:num_samples
        random_phases = 2*pi*rand(1, num_modes);
        random_amplitudes = rand(1, num_modes);
        random_amplitudes = sqrt(total_intensity) * random_amplitudes / norm(random_amplitudes);
        
        coeffs = random_amplitudes .* exp(1i*random_phases);
        coefficients(i, :) = coeffs;
        
        superposition = zeros(rows, cols);
        for j = 1:num_modes
            superposition = superposition + coeffs(j) * LG_modes(:, :, j);
        end
        superpositions(:, :, i) = abs(superposition).^2;
    end
end
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import eval_genlaguerre
#


def main():
    num = 64
    l_values = [0, 1]  # 修改变量名
    p_values = [0, 1]  # 修改变量名
    n = 1000

    threshold = 0.95

    # 生成LG光学数据
    lg_modes = lg_modes_func(num, l_values, p_values)
    datas = generate_random_superpositions(lg_modes, n)

    # PCA和绘图
    pc_exp, components_exp, scores_exp, n_min = pca(datas, threshold)
    plot_results(pc_exp, components_exp, scores_exp, n_min, num, datas.shape[0])

    print('\nCumulative Explained Variance:')
    for i in range(n_min):
        print(f'PC{i + 1}: {components_exp[i]:.1f}%')
    print(f'First {n_min} PCs cumulative: {sum(components_exp[:n_min]):.1f}%')


def pca(data, threshold):
    rows, cols, num_samples = data.shape

    # 重塑数据矩阵 - 确保与MATLAB相同的维度
    data_matrix = data.reshape(rows * cols, num_samples).T

    # 中心化数据 - 使用实数部分
    mean_data = np.mean(data_matrix, axis=0)
    data_centered = data_matrix - mean_data

    # 计算协方差矩阵 - 使用实数部分
    covariance_matrix = np.real(data_centered.T @ data_centered) / (num_samples - 1)

    # 特征值分解 - 确保使用实对称矩阵
    eigenvalues, eigenvectors = np.linalg.eigh(covariance_matrix)

    # 按特征值降序排序 (eigh返回的是升序)
    eigenvalues = eigenvalues[::-1]
    eigenvectors = eigenvectors[:, ::-1]

    # 计算累计解释方差
    total_variance = np.sum(eigenvalues)
    cumulative_variance = np.cumsum(eigenvalues) / total_variance

    # 找到达到方差阈值所需的主成分数量
    n_min = np.where(cumulative_variance >= threshold)[0]
    n_min = n_min[0] + 1 if len(n_min) > 0 else len(eigenvalues)

    # 提取所需的主成分
    pc = eigenvectors[:, :n_min]
    scores = data_centered @ pc
    components = 100 * eigenvalues[:n_min] / total_variance

    return pc, components, np.real(scores), n_min


def plot_results(pc, components, scores, n_min, n_size, data_rows):
    # 创建图形布局
    plt.figure(figsize=(15, 10))

    # 主成分图像
    n_plots = min(n_min, 6)  # 最多显示6个主成分
    for i in range(n_plots):
        plt.subplot(2, 3, i + 1)
        # 重塑并显示主成分
        pc_img = pc[:, i].reshape(n_size, n_size)
        plt.imshow(np.real(pc_img), cmap='jet')
        plt.axis('off')
        plt.title(f'PC{i + 1} ({components[i]:.1f}%)')
        plt.colorbar()

    # 如果主成分少于6个，绘制散点图在最后一个位置
    if n_plots < 6:
        plt.subplot(2, 3, n_plots + 1)
    else:
        plt.figure(figsize=(12, 8))
        plt.subplot(2, 2, 1)

    # 散点图 - 使用前两个主成分
    plt.scatter(scores[:, 0], scores[:, 1], s=20, alpha=0.6,
                color=[0.2, 0.4, 0.8])
    plt.xlabel('PC1')
    plt.ylabel('PC2')
    plt.grid(True)
    plt.axis('equal')
    plt.title('PC1 vs PC2')

    # 方差条形图
    if n_plots < 6:
        plt.subplot(2, 3, 5)
    else:
        plt.subplot(2, 2, 2)

    n_bars = min(10, n_min)
    plt.bar(range(1, n_bars + 1), components[:n_bars])
    plt.xlabel('Principal Components')  # 修改为英文
    plt.ylabel('Explained Variance (%)')  # 修改为英文
    plt.title('Variance per Principal Component')  # 修改为英文

    # 累计方差图
    if n_plots < 6:
        plt.subplot(2, 3, 6)
    else:
        plt.subplot(2, 2, 3)

    cumulative = np.cumsum(components)
    plt.plot(range(1, n_min + 1), cumulative, 'o-', linewidth=2)
    plt.axhline(y=95, color='r', linestyle='--', linewidth=1)
    plt.xlabel('Number of Principal Components')  # 修改为英文
    plt.ylabel('Cumulative Explained Variance (%)')  # 修改为英文
    plt.title('Cumulative Explained Variance')  # 修改为英文
    plt.ylim(0, 100)

    plt.tight_layout()
    plt.show()


def lg_modes_func(n_size, l_values, p_values):
    x = np.linspace(-3, 3, n_size)
    x_grid, y_grid = np.meshgrid(x, x)  # 改为小写
    r = np.sqrt(x_grid ** 2 + y_grid ** 2)  # 改为小写
    phi = np.arctan2(y_grid, x_grid)  # 改为小写

    w0 = 1
    lg_modes = []

    for l_val in l_values:
        for p_val in p_values:
            # 使用scipy的拉盖尔多项式
            ll = eval_genlaguerre(p_val, abs(l_val), 2 * r ** 2 / w0 ** 2)
            amplitude = (np.sqrt(2) * r / w0) ** abs(l_val) * ll * np.exp(-r ** 2 / w0 ** 2)
            phase = np.exp(1j * l_val * phi)
            mode = amplitude * phase
            # 归一化 - 使用与MATLAB相同的方法
            mode = mode / np.sqrt(np.sum(np.abs(mode) ** 2))
            lg_modes.append(mode)

    return np.stack(lg_modes, axis=2)

def generate_random_superpositions(lg_modes, num_samples):
    rows, cols, num_modes = lg_modes.shape
    superpositions = np.zeros((rows, cols, num_samples))
    total_intensity = 1

    for i in range(num_samples):
        random_phases = 2 * np.pi * np.random.rand(num_modes)
        random_amplitudes = np.random.rand(num_modes)
        # 归一化振幅 - 确保与MATLAB相同
        random_amplitudes = np.sqrt(total_intensity) * random_amplitudes / np.linalg.norm(random_amplitudes)

        coeffs = random_amplitudes * np.exp(1j * random_phases)

        superposition = np.zeros((rows, cols), dtype=complex)
        for j in range(num_modes):
            superposition += coeffs[j] * lg_modes[:, :, j]

        # 取模的平方得到强度
        superpositions[:, :, i] = np.abs(superposition) ** 2

    return superpositions


if __name__ == "__main__":
    main()